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Abstract 

(N 

' A grazing bifurcation corresponds to the collision of a periodic orbit with a switching 

manifold in a piecewise-smooth ODE system and often generates complicated dynamics. 
The lowest order terms of the induced Poincare map expanded about a regular grazing 
CN I bifurcation constitute a Nordmark map. In this paper we study a normal form of the 

■ Nordmark map in two dimensions with additive Gaussian noise of amplitude, e. We show 

that this particular noise formulation arises in a general setting and consider a harmoni- 
cally forced linear oscillator subject to compliant impacts to illustrate the accuracy of the 
map. Numerically computed invariant densities of the stochastic Nordmark map can take 
highly irregular forms, or, if there exists an attracting period-n solution when e = 0, be 
well approximated by the sum of n Gaussian densities centred about each point of the de- 
terministic solution, and scaled by -, for sufficiently small e > 0. We explain the irregular 
forms and calculate the covariance matrices associated with the Gaussian approximations 
in terms of the parameters of the map. Close to the grazing bifurcation the size of the 
invariant density may be proportional to y/e, as a consequence of a square-root singularity 
in the map. Sequences of transitions between different dynamical regimes that occur as the 
primary bifurcation parameter is varied have not been described previously. 



1 Introduction 

Deterministic, piecewise-smooth systems are commonly used as mathematical models of vibro- 
impacting systems [H EJ EJ H]. Specific examples include atomic force microscopy [5j E], gear 
assemblies [8], metal cutters [9j [10], and vibrating heat-exchanger tubes [HI [12]. Piecewise- 
smooth systems are characterized by the presence of co dimension-one regions of phase space, 
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termed switching manifolds, across which the functional form of the system changes. In the 
context of vibro-impacting systems, switching manifolds correspond to an impact or loss of con- 
tact. Motion that involves recurring impacts is often born in the collision of an attracting 
non-impacting periodic orbit with a switching manifold. This collision is known as a grazing 
bifurcation, Fig. [TJ Grazing bifurcations have been intensely analyzed and may be the cause of 
complex dynamics including chaos, see [13] and references within. In this paper we study the 
effect of noise on a common class of grazing bifurcations by analyzing a stochastic version of the 
induced return map. 

The effect of small noise on a well-understood deterministic system is often intuitive and 
qualitatively predictable but in certain situations may produce interesting new dynamics. For 
instance there are various mechanisms, such as coherence resonance, by which the addition of 
noise to quiescent systems induces relatively regular oscillations [HI [151 EE], as we ^ as more 
complicated dynamics such as mixed-mode oscillations [17]. Alternatively noise may suppress 
chaos [T3]. Studies of stochastic versions of some prototypical models are found in [TjJJ [2UJ l2~Tj 

123 [23]. 

In [2H [22], Griffin and Hogan study a one-dimensional piecewise-linear map with additive 
noise. For uniformly distributed noise they calculate basins of attraction and widths of invariant 
densities. With Gaussian noise the map resists such a precise analysis and instead numerical 
results are presented revealing that the two types of noise produce qualitatively similar behaviour. 
In a period-incrementing scenario, the invariant density of the stochastic map may undergo several 
transitions (as distinguished by the number of peaks in the density) as the noise amplitude is 
varied [26] . Numerical simulations of a two-dimensional ODE system with additive noise designed 
such that the relevant Poincare map is piecewise-linear, give essentially the same dynamics as the 
map. In [27J, numerical simulations of a noisy piecewise-smooth map are shown to exhibit many 
similarities to the output of a DC/DC converter. Noise- induced stabilization in one-dimensional 
discontinuous maps is described in [25] . 

Impact oscillators subject to noise have been studied from the point of view of characterizing 
small- amplitude, noise-induced oscillations. The early work of Dimentberg and Menyailov [29] 
concerns an unforced, damped, linear oscillator experiencing noise and instantaneous impacts 



/j, < /i = ^ > 




Figure 1: Phase space schematics illustrating a grazing bifurcation. The vertical line, E, 
represents a switching manifold. The bifurcation occurs when a periodic orbit intersects the 
switching manifold as a parameter, /i, is varied. Dynamics shortly after the bifurcation may be 
extremely complex and depend on the stability multipliers of the non-impacting periodic orbit 
and the precise nature of the grazing scenario [T3] . 
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with either one wall or two symmetrically placed walls. Zhuravlev's transformation [30| 0] con- 
verts the stochastic equation of motion with an impact rule into a single discontinuous stochastic 
differential equation. To then determine the invariant density of the position and velocity of the 
oscillator, in the case that impacts are slightly inelastic the system may be treated as a perturba- 
tion from a Hamiltonian system. Here stochastic averaging [5TJ E2J E5J leads to a simple stochastic 
differential equation for the energy. The work [29J has been extended to compliant impacts [51] , 
Hertzian contacts (35], and nonlinear oscillators [36J. The presence of forcing adds significant 
difficulty to the problem but has recently been investigated numerically and via stochastic av- 
eraging [53 EZ] and analyzed with series expansions and mean Poincare maps in the context of 
gear rattling [5B1 155] . 

Grazing bifurcations are typically studied by deriving and analyzing a Poincare map that 
is valid in a neighbourhood of the bifurcation. Different grazing scenarios yield fundamentally 
different Poincare maps. In two dimensions, a normal form of the Nordmark map is: 



x 



A 
A 



X 




" " 






y 


+ 


1 




) 




X 









y - 


- xVx _ 


+ 


1 



x < 
/x , x > 



xv, 



where 



.4 



r 1 
-6 



X = ±l 



;i-2) 



and t,5 Gl. This map applies to both regular grazing in three-dimensional, piecewise-smooth 
systems [40J, and grazing in vibro-impacting systems for which impacts are modelled as instan- 
taneous events with velocity reversal and energy loss [HJ H2]. The grazing bifurcation occurs at 
the origin when \i = 0. Since f 1 1.11) is a truncated series expansion centred about the bifurcation, 
(11.11) is primarily of interest for small values of x, y and \i. The construction of (11. ip is described 
in §2.11 Here we note that the fixed point of the left half-map of (11.11) is given by 
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and corresponds to the non-impacting periodic orbit. Assuming this periodic orbit is attracting, 
the fixed point, (11.31) . is admissible (i.e. x*^ < 0) when fi < 0. 

Here we describe the basic results and summarize the remainder of the paper. We analyze 
(11.11) with small-amplitude, additive Gaussian noise: 
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that is, £ is a two-dimensional Gaussian random variable with zero mean and covariance matrix 
0. In §2.21 we demonstrate that this noise formulation arises when small-amplitude Brownian 
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motion is added to the underlying differential equations and obtain an explicit expression for 
G. For a model of a linear oscillator experiencing impacts with a compliant support, in §2.31 we 
show that (11. 4p exhibits quantitatively the same dynamics as the stochastic differential equation 
model. Calculations for this section are deferred to Appendix [A3 

Figs. |2] and [3] show bifurcation diagrams of f ll.4p for two different choices of the parameter 
values. These were computed numerically by iterating (ll.4p over a range of fixed values of \i. 
As described in [25] , the noise blurs the underlying deterministic bifurcation diagram and causes 
impacts to occasionally occur when [i < 0. In §3.11 we analytically compute periodic solutions 
of (11. 4p with e = 0. These are indicated by blue curves in Figs. [2] and [3J For the parameter 
values of Fig. [2J in the absence of noise the maps exhibits a period-incrementing cascade, [4"3t IT3]. 
For every n > 1, there exists an interval of positive /j,- values for which period- n solutions are 
admissible and attracting (periods 3, 4 and 5 are indicated in Fig. [2j. These intervals are ordered 
by n, disjoint, and separated by small chaotic bands. As n — > oo, the intervals approach /i = 
from the right. 

When e > 0, dynamics may be roughly periodic (e.g. when /i = 0.01 in Fig. [2]), however, 
the noise smothers high-period solutions. As with the logistic map [SI 05], the highest period 
that is in some sense distinguishable, decreases with an increase in the noise amplitude. Figs. [2] 



x 




-0.12 1 1 1 1 1 1 ' 

-0.005 0.005 0.01 0.015 0.02 0.025 

/* 

Figure 2: A bifurcation diagram of the stochastic Nordmark map, (11.41) . with (r, S, x) = 
(0.5,0.05,1), e = 0.00025, and = I (the 2x2 identity matrix). The dots are numerically 
computed iterates of (II. 4p with transients omitted. The blue curves denote deterministic, admis- 
sible, attracting periodic solutions, computed analytically in §3.11 The red curves approximate 
one standard deviation of the invariant density from the blue curves, as obtained via linearization, 
see §EJ2J 
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and [3] suggest the existence of a stable invariant density in a neighbourhood of the underlying 
attractor, for any /j,. Throughout the paper we study this density but do not formally prove its 
existence or uniqueness. The invariant density may be approximately Gaussian about each point 
of a deterministic periodic solution of relatively low period. In §3.21 we derive the covariance 
matrices of these Gaussians. (The Gaussian approximation to the unimodal invariant density 
centred about the fixed point, (11.31) . for /i < 0, may be computed in the same way.) In §3.31 
we derive an approximation to these covariance matrices and in §3.41 investigate the effect of 
different O, (II. 5p . The near-Gaussian densities are discussed further in §4.11 The accuracy of the 
Gaussian approximation decreases as e is increased, §4.21 For relatively large e, the density is 
strongly non-Gaussian and may exhibit some interesting novel characteristics, §4.31 Finally, we 
note that when e = attracting solutions may coexist. For instance in Fig. [3l for small /i < 
the fixed point (I1.3P coexists with a period-3 solution, and for a relatively wide range /i- values in 
the middle of the figure, this period-3 solution coexists with either a period-2 solution or what 
we numerically observe to be a chaotic attractor. When e > 0, numerics indicate the existence 
of a single invariant density with peaks at the deterministic attractors, §4.41 Finally conclusions 
are presented in §5J 




Figure 3: A bifurcation diagram of (11.41) . with (r, 5, x) — (0-2, 0.5, 1), e = 0.01, and = 1. The 
meaning of various figure elements are explained in the caption of Fig. [2j 
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2 Derivation of the stochastic Nordmark map 



The purpose of this section is to justify the manner by which we have incorporated noise into 
the Nordmark map, (11.41) with (jl.5p . and derive an explicit expression for in terms of an 
underlying stochastic differential equation. Certainly noise may enter a vibro-impacting system 
in a variety of different ways. Noise in external forcing may be represented by a deterministic 
signal plus a weak random component modelled by white or coloured noise. Model parameters 
may be treated stochastically. Also impact events may be modelled randomly (indeed experiments 
reveal impacts may induce high frequency oscillations that quickly decay and vanish [IB])- The 
manner by which randomness is modelled affects the nature of the induced Poincare map. Below 
we consider additive, white Gaussian noise and derive (II .4p . 



2.1 Construction of the deterministic Nordmark map 

To describe the construction of the Nordmark map it is necessary to introduce a general system 
that exhibits a regular grazing bifurcation. To achieve this in a general setting requires consider- 
able space in order to carefully state assumptions, non- degeneracy conditions and explain their 
meaning. For brevity we omit most of these statements - details are given in [T3j . 

Consider an ODE system in M 3 , u = (u,v,w), with a single smooth switching manifold, 

£ = {u|w = 0}, (2.1) 

that we write as 

f /W(u;iy), u<0 

where and are smooth functions extendable beyond their respective half-spaces, and 
T] is a parameter. Suppose when r/ = there exists an attracting periodic orbit that intersects 
u = but is otherwise contained in {u | u < 0}. To negate the possibility of sliding motion [T3j . 
and specify a direction of flow, suppose 

sgn (eJfM (0,v,w;ri)) = sgn (ejf™ (0, v, w; r/)) = sgn{v) , (2.3) 

in a neighbourhood of u = 0. 

We introduce two Poincare sections: II is a generic cross-section transversal to the flow lying 
entirely in the left half-space, and 

n' = {u | v = 0} , (2.4) 

see Fig. HJ To analyze oscillatory dynamics one can study the induced return map on II, or a 
return map on II' using intersection points that for u > are virtual, as explained below. The 
Poincare map on II is useful for comparing with numerical simulations. However, we use the 
section II', because, although some iterates are virtual, the map takes a simpler form and is more 
amenable to analysis. 

To obtain the return map on II', consider the trajectory from a point, uo, on II, as governed 
by (12.21) . and suppose that this trajectory intersects £ at u 2 and u 3 before returning to II at some 
u 5 , as shown in Fig. HI We use the left half-flow, f^ L \ to follow the trajectory beyond u 2 into the 
right half-space until intersection with II' at u x . Similarly we use to follow the trajectory 
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backwards from 113 to U4. The map from Ui to U4 defined in this manner, call it D, is known as 
a discontinuity map and takes the form 
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ui < 
, ui > 



(2.5) 



for some scalar constant c. This concept was originally introduced by Nordmark [HJ H2] in 
regards to instantaneous impacts and is nowadays a fundamental tool for analyzing a variety of 
grazing phenomena [4T)l Wf\ |4~8| I49| 150] . The discontinuity map is nonsmooth and for u\ > 
corresponds to motion on an 0(-y/tti) time scale, where Ui is assumed to be small. In contrast, 
the induced Poincare map from 114 to Ug, as governed by the left half- flow, captures the global 
dynamics on an 0(1) time scale. This global map, which we call G, is smooth and thus may be 




Figure 4: A schematic for regular grazing in three dimensions. We show intersections of a 
trajectory of the system (12.21) with a generic Poincare section, II, the alternate Poincare section, 
IT, (12.41) . and the switching manifold, S, (12. ip . Since the trajectory intersects the switching 
manifold at u 2 , the points Ui, U4 6 IT' are found by solving the left half-system, f^ L \ in the right 
half-space. The discontinuity map is u 4 = D(ui), (I2.5p . The global map is u 6 = G(u 4 ), (12. 6p . 
The normal form Nordmark map, (11.11) . is given by G o D, after the coordinate change (12.71) and 
truncation. 
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written as 

+ bri + 0(\u 4 ,w i ,r ] \ 2 ) . (2.6) 

for some constant matrix A and vector b. The Nordmark map (II. ip is then obtained by: (i) 
forming the composition GoD, (ii) transforming to normal form via the linear change of variables, 
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and (iii) omitting higher order terms in x, y and /i. The Nordmark map ignores linear terms that 
are dominated by the square- root term. For a careful consideration of these terms see |51j . 
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2.2 Addition of noise 

The stochastic differential equation: 

du = { /^(u; t \ « > } dt + £B(U] V) dW{t) ' (2 ' 9) 

where W(t) is a vector of independent standard Brownian motions, B is matrix with smooth 
dependency on u and i], and < e < 1, represents the system (12. 2p in the presence of small- 
amplitude, white, Gaussian noise. In this section we extend the construction described in the 
previous section to accommodate the noise term and obtain (jl.4p . We assume e is small enough 
that the basic global structure of the system is not destroyed by the noise and sample paths of (12. 9p 
likely stay near to the associated deterministic solution, over an 0(1) time frame. Theoretical 
results on the existence and uniqueness of solutions to stochastic differential equations with 
discontinuous drift may be found in [52], and references within. 

Our strategy is to obtain stochastic versions of the global map, G, (12. 6p . and the discontinuity 
map, D, (12. 5p . and compose them to arrive at the stochastic Nordmark map, (II. 4p . We first 
discuss the stochastic global map and consider solely the smooth left half-system of (12. 9p . Issues 
relating to switching are discussed below. 

When a sample path of the left half-system of (12 .9p passes through II', if the noise is not 
purely tangent to II', with probability 1 the sample path has multiple intersections with II', 
albeit over an extremely short time frame. We define the stochastic global map by the first 
reintersection, u 6 , of the sample path with II' after one large excursion. (An alternative is to use 
the mean passage time [52]. For stochastic Poincare maps for systems with coloured noise, see 
[54].) Multiple intersections occur sufficiently quickly that this specification does not affect the 
lowest order noise term in the map [55] . 

For e > 0, the quantity U6 is stochastic and to calculate the stochastic global map means to 
compute the probability density function for u 6 . We can do this to 0(e) by applying standard 
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asymptotic results usually stated in the context of an exit problem. Here we follow the sample 
path methodology of Freidlin and Wentzell [55] , see also [56| 157] . 
For the general stochastic differential equation 

dv = /( v ) dt + eB(v) dW{t) , v(0) = v , (2.10) 

where / and B are C 2 , by Theorem 2.2 of Chapter 2 of [55] we can write 

v(t) = v(°\t) + ev (1 )(t) + o(e) , (2.11) 

where v( ) is the deterministic solution to (I2.10p with e — 0, and v^ 1 ) is the solution to the 
time-dependent Ornstein-Uhlenbeck process: 

rfv (1) = AJ(v (0) )v (1) + £(v (0) ) , v (1) (0) = . (2.12) 

We then have the following lemma. 



Lemma 1 (Theorem 2.3 of Chapter 2 of Freidlin and Wentzell 

Let V be a domain containing v and let be the first passage time o/v^ from T>. Suppose 
dT> is differentiate at v^(T^), let p denote the exterior unit normal vector to dT> at this point, 
let q = f (v(°)(T(°))), and assume p T q ^ 0. IfT denotes the first passage time o/v from V, then 

T = T^- pTvil) y i0)) e + o(e), (2.13) 
p' q 

v(T) = V (°)(T(°)) + (j - ^ v( 1 )(T(°))e + o(e) . (2.14) 

V p?/ 

The key quantity, v^(T^), is stochastic, and by ( 12. 12H . 

v« (t) = jf Exp ^ D v / (v (°) (r)) rfr^) B (v<°) (s)) . (2. 15) 

Consequently v^(t) is a Gaussian random variable with zero mean and covariance matrix, call 
it Q(t), given by 



where 



Q(t)= [ H{s,t)H{s,t) T ds , (2.16) 
Jo 

H(s,t) = Exp (f s D vf (v (0) (r)) d/) B (v^(s)) . (2.17) 

Therefore, to 0(e), the location of first exit, v(T) G <9P, is Gaussian with mean, v^(T^), and 
covariance matrix 

(V^W (o) )f/-^y . (2.18) 

V pqj V p?/ 

In the context of the global map, / = f^ L \ and &D represents IF. By Lemma [TJ 



w 6 



G(ut, w A - rj) + eCc + o(e) , (2.19) 
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where G is the deterministic global map, ( 12. 6p . and £g is a zero- mean, two-dimensional Gaussian 
random variable. Herep= [0, — 1, 0] T and q = [0, — 7l, Cl] T +0(|m4, W4, rj] 1 ), where ([0,0,0] T ;0) 
[0, — 7L)Cl] T - Consequently, by (12 . 1 8j) . the covariance matrix for £ G is 



^U 12 + W13 



~1L 



U)\2 + W13 



^T^22 + ^~ U 23 + <^33 



+ 0(|M 4 ,t«4,^| 1 ' 



(2.20) 



where Uy denotes the (i, j) -element of fi(T ( - ^), (12. 16ft . and denotes the period of the grazing 
periodic orbit at 77 = 0. 

A stochastic version of the discontinuity map is not as straightforward. For instance, consider 
the motion from U2 to 113, Fig. HI In the absence of noise, 113 is easily calculated using the right 
half- flow. However, with noise acting non-tangentially to S, with probability 1, a sample path 
from U2 to U3 spends a nonzero amount of time in the left half-space. Thus motion from U2 to 
u 3 is in part governed by the left half-system, albeit weakly. Moreover, to calculate the motion 
from Ui to u 2 and from u 3 to u 4 , we must solve the stochastic differential equation backwards 
in time. 

Here we argue a bound for the noise in the discontinuity map. Intuitively, the noise in the 
discontinuity map is greater for larger values of U\ > 0. For fixed e > 0, consider u\ 3> e, but recall 
we assume U\ is small. In this case the excursion from Ui to U4 is well-approximated by three exit 
problems for which Lemma [1] is applicable, because the weak noise does not significantly disrupt 
the three pieces of the excursion. The discontinuity map acts over an 0(- v /u 1 ~) time- frame, and 
note that a consequence of Lemma [1] is that for small t, the standard deviations of components 

of the exit point are 0{ey/t). Therefore we expect deviations in u 4 to be at most O (euj\ and 
write 



"1 

W4 



D^w^rj) + e£ D + o(e) 



(2.21) 



where £d < O yui 

When we compose (12.191) and (12.211) . assuming u\ is small, £g dominates £d because £g = 
0(1). As in §2.11 to this composition we apply the change of variables (12. 7p . and omit higher 
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order terms, producing ( II ,4p . Moreover, the noise term in ( II .4p is ^ 
and therefore the covariance matrix for £ is given by 
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2.3 A linear oscillator with compliant impacts 

As an example, we consider a harmonically forced linear oscillator experiencing compliant impacts 
with a massless, prestressed support. The system is depicted in Fig. [5] and we model it with the 
non-dimensionalized equation of motion, 



u 



-k osc (u + 1) - b osc ii , u<0 
-k osc (u + 1) - (b osc + b supp )u - A; supp (w + d) , u > 



Fcos(t) +ev(t) , (2.23) 
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where u(t) is a standard Gaussian white noise process. Here u denotes the position of the mass. 
A rigid stop placed at u = prevents the support moving left of this value and prestresses the 
support by a distance d > 0. We assume that while the mass is not in contact with the support, 
the support is at rest at u — 0. Whenever the mass arrives at u — 0, the mass and support then 
move together for u > until a return to u = at which point the stop catches the support 
and the mass returns to non-impacting motion. The setup is identical to a scenario given in [58] . 
Experiments involving compliant impacts are described in [591 EH] • We incorporate randomness 
in the model by adding white-noise to the forcing term. 

b 2 

To simplify our discussion we assume, < k osc < 1, so that, in particular, the system is 
sub-resonant. We treat the forcing amplitude, F, as the primary bifurcation parameter. When 
e — 0, for small F the system settles to a stable non-impacting periodic orbit of period 2n. The 
periodic orbit undergoes a grazing bifurcation at 

F graz = v^L + U-^osc) 2 . (2.24) 

A typical bifurcation diagram of f l2.23|) with e > is shown in Fig. [6]A. Fig. |6]-B is a bifurcation 
diagram of the stochastic Nordmark map, (11.41) . with parameter values corresponding to the 
impact oscillator for panel A. As detailed in Appendix |AJ the values r, 5 and x, are given by 
the deterministic results of §2.11 and is given by (12.221) . As expected, the two bifurcation 
diagrams exhibit the same qualitative structure. To quantitatively compare (I2.23P and ( II .4p . 
for each we have computed the corresponding bifurcation diagrams on II', panels C and D. As 
expected these bifurcation diagrams are practically indistinguishable, justifying our use of the 
map (11.4jl . A more precise comparison is beyond the scope of this paper. 

3 Approximating invariant densities by a linear analysis 

In this section we consider the case where the invariant density can be approximated by a 
combination of Gaussians and calculate the related covariance matrices. The calculations reveal 
instances of larger variability and situations where the Gaussian approximation may break down. 

We consider the case that (11.41) with e = has an attracting periodic solution, {(xj, yi)}i=o,...,n-i, 
with Xi 7^ 0, for all i. Then with small e > 0, numerical simulations indicate that (11.4j) has a 




u = —1 u = —d u = 



Figure 5: A schematic of the forced vibro-impacting system modelled by ( 12. 23ft . The equilibrium 
positions of the mass and support are u = — 1 and u = —d respectively. 
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stable invariant density with high peaks near each (xi, yi). For sufficiently small e > 0, the bulk of 
this invariant density is distant from the switching manifold. Consequently in this case the non- 
smooth effect of the switching manifold is not seen. About each (x,, yi) the invariant density has 
a roughly Gaussian shape that may be computed by linearizing the map about (xj, yj). Therefore 
the invariant density of (11.41) is well-approximated by the sum of n Gaussian densities centred at 
each (xi, yi), and scaled by - such that the density is normalized (see Fig. [7] for an example). In 
§3.11 we derive an appropriate n th iterate of (11. 4p and solve for (xi,yi) when e = 0. In §3.21 we 
calculate the covariance matrices of the Gaussians. Further calculations and approximations are 
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Figure 6: Panels A and C show numerically computed bifurcation diagrams of the impact 
oscillator model (I2.23H with (A; osc , 6 osc , /c supp , 6 supp , d) = (5,0.5,10,0,0.1) and e = 5 x 10~ 5 . In 
panel A, the values are intersections of long-time dynamics with II = {(u,ii) \ ii = 0, u « 
—2}. Panel C shows intersections with IT. For u\ > these intersections are virtual and were 
obtained by appropriately solving the non-impacting equation of motion for u > 0. Panel B is a 
bifurcation diagram of the stochastic Nordmark map, (II. 4p . using parameter values corresponding 
to the system in panel A as calculated by the formulas given in the text. Specifically, (r, 5) ~ 
(0.07264,0.04321), x = 1, and (6> n , 12 , 22 ) ~ (662.6,-7.450,28.29), where % denotes the 
element of 0. Panel D shows the points of panel B under the inverse of the coordinate change 

(EZD- 
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described in §3.31 and §3.41 



3.1 Maximal periodic solutions 

Stability of a periodic solution, {{x ir yi)}i=o,...,n-ij of (11.41) with e = may be determined in 
a standard manner from the multipliers of the matrix formed from the product of Jacobians 
evaluated at each (xj, y^). For each point of the solution for which Xi > 0, the product accumulates 

a term of order x i 2 , which is strongly destabilizing since we assume the Xi are small. For this 
reason we only consider periodic solutions for which exactly one iterate satisfies Xi > 0. We 
expect no other periodic solutions to have a significant effect on long-term dynamics near the 
grazing bifurcation. Indeed this is consistent with numerical results. This observation has been 
noted previously and such solutions are known as maximal periodic solutions [131 [61] . 

Suppose that the next n—1 iterates of (11.41) from a point (xq, yo) with xq > satisfy X{ < 0. 
Then the n th -iterate is given by 



Xn 


= A n 


x 






_ yo - Xa/^o . 



n-l 



+ (/ + A + ■ ■ ■ + A n ~ l ) 
+ ■■■ + A n ~%) , 




1 



(3.1) 



where the n random vectors, £j, are independent. The covariance matrix of the sum of n inde- 
pendent Gaussian random variables the sum of the covariance matrices of each, therefore 



+ au-2 + ■■■ + A n -% ~ e (n) , 

where ^ n - ) is a zero-mean Gaussian random variable with covariance matrix 



n-l 



(3.2) 



(3.3) 



i=0 



Hence we can write (13.11) as 



Xn 


= A n 


x 






_ yo - xV^o _ 



+ 6 (n V + ei 



where 



A 



+ A 



n-l> 




1 



b {n) = (/ 

Fixed points, (x* {n) , y* {n) ), of fl33D with e = 0, satisfy 

a 22 b 1 ) X Vx^) + ((1 - a u )b 2 + a 21 b 1 ) x*^ 



y 



in) 



,a 12 b 2 



a 22 )bt + a 12 b 2 



a 2 xVx*^> + ((1 - on)(l - a 2 2) - ai 2 a 2 i) x* 



(n) 



where 



A r ' 



a n 0,12 
a 2 \ a 22 



1 - 022)^1 + ai 2 b 2 
, &<»> = 



t>2 



(3.4) 
(3.5) 

(3.6) 
(3.7) 

(3.8) 
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The stability multipliers associated with (x< n \y*^) are the eigenvalues of the matrix 



1 


" 




an - 


0-12X 


Ol2 




hi hi 




lVx*( n ) 




-x 


1 




a>2i ~ 




«22 




hi hi 


2v / £*( n ) 











(3.9) 



In summary, the point of a maximal periodic solution with x > 0, (x 



*(n) ^ *(n) 



is given im- 
plicitly in terms of \i by (13. 6p and (13.71) . The remaining points are easily obtained by iterating 
(11.41) from this point. The periodic solution is stable when both eigenvalues of K(x*^), (13. 9p . 
lie inside the unit circle. The blue curves in Figs. [2] and [3] denote the x- values of stable, admis- 
sible, maximal periodic solutions (up to period 5) and were computed by evaluating the above 
expressions. 



3.2 Calculation of Gaussian invariant densities 

Here we determine the covariance matrices for the Gaussian approximations. We first compute 
the covariance matrix for the Gaussian centred at (x*^ n \ y*^), then give an iterative formula for 
the remaining covariance matrices. 



The linearization of the n th -iterate, (13.41) . about the fixed point, {x' 



IS 



K Zj 



(n) 



where K = K{x< n \^)), lEnty . and 



(3.10) 



X - £* (ri) (/i) 

y-y< n Hfj) 



(3.11) 



denotes the displacement from the fixed point. Since (13.101) is linear and has a Gaussian 
distribution (with zero mean and covariance matrix 6^, (E3D), fEUOl) maps one Gaussian prob- 
ability density function of z- values, into another [621 [63] . Specifically, if the j th density under 
(I3.10p is Gaussian with zero mean and covariance matrix s 2 Aj, then 



A i+ i = KAjK T + 6 (n) 



(3.12) 



If both eigenvalues of K lie inside the unit circle, then the fixed point (x*^ n \ y*^>) of (13.41) 
is attracting and (I3.10p converges to a unique, zero-mean, Gaussian, invariant density with 
covariance matrix, e 2 A, satisfying 



A = KAK T + 9 (ri) 



(3.13) 



Equation (13.131) is a Lyapunov equation that cannot be solved for A by elementary matrix 
algebra operations [22]. However we can break (I3.13P into its constituent scalar equations. We 
write (13.131) as 



-^11 ^12 

^12 A22 



hi hi 
hi hi 



-^11 ^12 
A12 A22 



hi hi 
hi hi 



+ 



"11 &11 
n{n) n(n) 
°11 °11 



(3.14) 
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which may be broken up and rewritten as 







' A n " 






Al2 


= M 


A12 


+ 


^12 


A22 




A22 




/)(«) 
. ^22 



(3.15) 



where 



Then 



M 



k 2 

ft 12 



Aft 2fciifci2 

fcllAfcl fcllfcjffl + k 12 k 21 k 12 k 22 
2k 2 \k 22 



k 2 

ft 21 



k 2 

^22 



An 
A12 
A22 



[I-M)- 



0$ 
0% 



(3.16) 



(3.17) 



A direct calculation reveals 
det(J-M) = (det(iT) 



trace(iT) + l) ( det(K) + trace(K) + l) (l - det(K)) , (3.18) 



thus det(J — M) = if and only if K has an eigenvalue on the unit circle. Moreover, if under 
parameter change, an eigenvalue of K approaches the unit circle, A tends to infinity. 

A, given by explicitly by (I3.17P with (13.1 6p and (13. 9p . is the covariance matrix for the Gaussian 
approximation to the invariant density about {x*^ n \ y*^). The covariance matrix relating to the 



i th iterate of 



{x* {n) ,y* {n) ] 



call it AM may be calculated via the identity 

A (i+i) = aa m a t + ) 

(0 



(3.19) 



(where A^ n ^ = A^ = A), which is valid for % — 1, . . . , n — 1. Let A n denote the (1, l)-element of 



A«. The red curves in Figs. |2] and |3] denote a distance eX[{ from the i th iterates of the periodic 
solution (the blue curves) and correspond to one standard deviation from the mean. The AS? 
may differ significantly for different i as in Fig. |6j 



3.3 An approximation to the covariance matrices 



Recall, the stability multipliers of a maximal periodic solution with a point, (x* {n \y* {n) ), in the 
right half-plane are the eigenvalues of K, (I3.9p . Near grazing, x*^ is small, thus if this periodic 
solution is attracting, the elements, a\ 2 and a 22 , of A n , must also be small. In view of (13.91) . it is 

reasonable to suppose |ajj| < O (V x*( n ^ , for each a^. Then from (13. 16 j) . 



and thus 



M 



"12 

4 x *(n) 
«12«22 

a 22 



(I- M)- 1 = 1 + 



4 x *(n) 










1 - 



+ O (Vx*( n )) 



7 2 
-'12 



a \2 

ai2d22 

a 22 
4a ..(n) 










(3.20) 



+ O (Vx*wj 



(3.21) 
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Then by (l3~T7j) . 



a = e (n) + 



(n) 



'12 



a 



12 a 12«22 



ai 2 a 22 a 



22 



(3.22) 



Equation ( 13 .22 j) indicates the deviation of A (the covariance matrix for the Gaussian invariant 
density about (x*^ n \y*^)) from 0( n ), (13. 3p . For large n we often have | a^- 1 Vr^, in which 



case A ps 0( n ). Alternatively, if the periodic solution is weakly attracting, then by ( 13. 181) . the 



denominator, Ax*^ 



12' 



is small. By (13. 22 p . in this case A deviates significantly from Q^ n \ as 



seen in Figs. [2] and |3] where the red curves are relatively distant from the associated blue curves. 



3.4 Effects of different 9 

Both eigenvalues of A lie inside the unit circle (assuming the grazing periodic orbit is attracting) 
and therefore the infinite series, 



Q (oo) =i^ A iQ ( A ry ^ 



(3.23) 



i=0 



converges. Thus for large n, by (13. 3p . 0(°°) may provide a good approximation to 0( n ). Here we 
derive a simple expression for 0* 00 ). Equation (I3.23P is equivalent to the Lyapunov equation, 



= 9 (oo) - AB {co) A T . 
By expanding the individual elements of this equation and solving we obtain 



(3.24) 



°11 


1 


^12 
fl(°°) 
. ^22 


~ A 



1 + 6 
-t6 

8 2 + 5 3 



2r 

1 - r 2 - 5 2 
2t5 2 1 



1 + 6 
-tS 

8-r 2 + r 2 5 





' 011 " 








6*22 



where 



A 



(3.25) 



(3.26) 



(5 - r + 1)(5 + r + 1)(1 - 5) > . 

We can use (13.251) to investigate the effects of different 0. In order to keep the noise amplitude 
controlled by e, here we fix trace(0) = 2. Substituting 6>n = 2 — 8 22 into (13.251) yields 









U 12 




1 




A 


U 22 







1 (1 



25 2 



8 

6 + t8 12 



1 



t9 12 ) 

2 



r~ - 8 2 )9 12 ] 



2 ) 
622 



(3.27) 



where Q\ 2 and Q 22 are constrained by 

e\ 2 + (e 22 - 1) 2 < 1 



< e 22 < 2 . 



(3.28) 



From (I3.27P we learn that the marginal variance, 0^° , depends only on r, 5 and the correlation 
of 0. Furthermore, if A is not too small, then if 5 is small (as in Fig. ED, e { 22 ] « 22 , and if r is 
small (as in Fig. [3]), 0^ may assume only a small range of values. 
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Figure 7: The invariant density of the stochastic Nordmark map, (jl.4p . using the same parameter 
values as Fig. [2] with also fi = 0.005. Here the noise amplitude is sufficiently small (e = 0.00025) 
that the density is approximately the sum of four Gaussians, scaled by \. The value of the 
density is indicated by color (dark red - the maximum value of the density; dark blue - zero). 
The vertical gray line denotes the switching manifold, x — 0. 

4 Invariant densities 
4.1 Near Gaussian densities 

As described at the beginning of §31 if (ll.4p with e = has a periodic solution, {(xj, yi)}i=o,...,n-i, 
satisfying Xj ^ for each i, and no other attractor, then for small e > we observe an invariant 
density that is well-approximated by the sum of n Gaussians centred about each (xi,yi), and 
scaled by -. The size of each Gaussian is 0(e). The precise shape of each Gaussian approximation 
is given by the results of §3.21 Specifically its mean is (xi,yi), and its covariance matrix, , is 
given iteratively by (I3.19p . where A^ ^ is given by (13.171) . 

An example is shown in Fig. [71 This figure, and others below, shows the invariant probability 
density function where the value of the function is indicated by colour (dark red - the maximum 
value of the density; dark blue - zero). Following [25], we assumed the system is ergodic and 
approximated invariant densities by a two-dimensional histogram computed from 10 8 consecutive 
iterates of a single orbit. 

With the parameter values of Fig. [7J iterates of (11.41) follow close to the underlying determin- 
istic period-4 solution. Iterates step from the neighbourhood of each point in the order left to 
right. While x < 0, iterates approach the deterministic fixed point of the left half-map, (II. 3p . 
which for \i > lies in the right half-plane and is virtual. Near the right-most point of the 
periodic solution, for which x > 0, the square-root term takes effect and the next iterate lies 
near the left-most point. The left-most Gaussian has the greatest spread because it corresponds 
to iterates that have just experienced the strongly expanding square-root term. Moving left to 
right, the size of each Gaussian decreases because the left half-map of f ll.4j) is contracting, and 
indeed we have the explicit formula (I3.19P . Note that since the left-most Gaussian is most widely 
spread, the value of the density at its peak is substantially less than at the other three peaks. 
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fraction of iterates with x > 
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0.8 



""0 0.0005 0.001 0.0015 0.002 

e 

Figure 8: The fraction of points of an orbit (with transients decayed) of (II. Ah with x > that 
return to the right half-plane in exactly four iterations of (\lAh . The parameter values are the 
same as in Fig. [2j Each data point was computed numerically from 10 7 iterations of one orbit. 
The invariant density when e = 0.00025 is shown in Fig. [7J and when e = 0.00075 in Fig. |9j 



4.2 Weakly non- Gaussian densities 

As the value of e is increased, the above approximation to the invariant density worsens because 
with increasing probability iterates fall on the "wrong" side of the switching manifold. For an 
arbitrary orbit, let Xj denote the number of iterations that the j th -point of the orbit satisfying 
x > takes to return to the right half-plane. Fig. [8] is a plot of the fraction of Ij that are 
equal to n, where n is the period of the deterministic periodic solution, against e, using the same 
parameter values as Fig. [21 and /j, = 0.005. 

A key feature of Fig. [8] is that, roughly speaking, this fraction begins to descend rapidly over 
a narrow interval of values of e (around e = 0.0006). Two main factors influence the location of 
this critical interval. If the deterministic periodic solution is weakly attracting, then the matrix, 
K, (I3.9p . has an eigenvalue near the unit circle, and thus by (I3.18P the matrix, / — M, is almost 
singular. Therefore by (13.171) the invariant density is broad and so the values of e in the critical 
interval are relatively small. Second, the values of e in the interval are roughly proportional to the 
distance of the deterministic periodic solution from the switching manifold. These observations 
explain the breakdown of the Gaussian approximation visible in Figs. [2]and[3j Near a saddle- node 
or period-doubling bifurcation the strength of attraction of the periodic solution is weak and the 
Gaussian fit is poor. Near a border-collision bifurcation the periodic solution is in close proximity 
to the switching manifold and again iterations often stray from the Gaussian approximation. 

When e = 0.00025, as in Fig. [7J iterations follow the deterministic period-4 pattern almost 
exclusively and the Gaussian approximation is a good fit to the numerically computed invariant 
density, as described above. When e = 0.00075, as in Fig. [HI iterates follow the period-4 pattern 
about 96% of the time. In this case the invariant density, shown in Fig. [9], exhibits four clear 
peaks near the deterministic solution, but level curves are relatively non-elliptical indicating that 
the density is not as well approximated by Gaussians. A relatively high fraction of iterates lie 
very near the switching manifold. Due to the square-root singularity, (11.41) maps a small region 
just to the right of the switching manifold to a large region that including points relatively distant 
from the deterministic periodic solution, say around (x,y) ~ (—0.02,0.005). Consequently the 
invariant density displays a nonlinear character. 
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Figure 9: The invariant density of fll.4j) when /x = 0.005, e = 0.00075 and the remaining 
parameter values are as in Fig. |2j 




x 



Figure 10: The invariant density of ( 11. 4ft when \i = and the remaining parameter values are 
as in Fig. [2j The dashed line is the slow manifold for the fixed point of the left half-map of ( II ,4p . 
Since the left half-map is linear, this manifold coincides with the eigenspace for the stability 
multiplier of greatest magnitude. 

4.3 Strongly non-Gaussian densities 

If e is large relative to the distance of the deterministic attractor from the switching manifold, or 
if this distance is zero, the invariant density takes a strongly non-Gaussian form. As an example, 
Fig. d0] shows an invariant density of ( 11.41) at the grazing bifurcation (/i = 0). Since we have a 
thorough understanding of the behaviour of map when e = 0, the shape of the density can be 
explained heuristically. Due to the square-root singularity, the map (11.41) flings points with x > 
(and y < X\fx — rx) relatively far into the left half-plane. While x < 0, subsequent iterates are 
governed by the left half-map and slowly step towards the fixed point of this map, which, when 
\x = 0, is located at the origin. For Fig. [101 the stability multipliers associated with the fixed 
point are real and distinct (approximately 0.14 and 0.36). Thus iterates first step quickly to the 
eigenspace corresponding to the multiplier of greater magnitude (0.36) (indicated in Fig. [TO]) , 
then head into the origin. This explains the angle at which the invariant density protrudes from 
x — 0. Since noise is present, iterates eventually pass over the switching manifold, and the 
process repeats. Since the noise amplitude is 0(e), iterates with x > are O(e). Due to the 
square-root singularity, the size of the invariant density is therefore 0(y/e), Fig. [TTJ In addition, 
when /i = the average number of iterations taken for points of an orbit to return to the right 
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Figure 11: Standard deviations of the marginal densities of the invariant density (shown in 
Fig. Unifor e = 0.00025) of (11.41) for different values of e when /i = and the remaining parameter 
values are as in Fig. |2l Each data point was computed numerically from 10 7 iterations of one 
orbit with transients decayed. 




x 

Figure 12: The invariant density of f ll.4jl when \l = 0.15 and the remaining parameter values 
are as in Fig. [3] 



half-space decreases as e is increased. 

4.4 Densities resulting from coexisting attractors 

If, in the absence of noise fll.4j) has coexisting attractors, with noise we expect there to be a single 
invariant density involving peaks near both deterministic attractors. Such a density is shown in 
Fig. H2J The relative weighting of the density about each deterministic attractor depends on 
the relative strengths of attraction and the relative sizes of the basins of attraction. Orbits may 
remain in the neighbourhood of one attractor for a long period of time before sufficiently extreme 
noise pushes the orbit into the neighbourhood of a different attractor. Similar behaviour has been 
described for smooth maps [54"] . 
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5 Conclusions 



In this paper we studied the Nordmark map with additive, Gaussian noise, (II .4p . by analyzing 
two-dimensional invariant densities that were computed numerically and in certain situations 
approximated by a sum of scaled Gaussian probability density functions. 

In £j2]we derived (11. 4p from a general three-dimensional stochastic differential equation, ( 12.91) . 
that, in the absence of noise, exhibits a regular grazing bifurcation. The map (11 .4p is the compo- 
sition of a smooth map encapsulating global dynamics, (the global map, (12. 6p ) and a piecewise- 
smooth local map (the discontinuity map, (I2.5P ) that represents the correction to the global 
map induced by the dynamics for u > 0. The noise term in ( 12. 6 p was calculated using standard 
stochastic asymptotics theory, §2.21 Near grazing the discontinuity map corresponds to dynamics 
occurring over a much shorter time period than the global map and for this reason the effect 
of randomness in the discontinuity map may be neglected. Moreover, with this assumption we 
derived the covariance matrix, O, for the noise in (jl.4p explicitly in terms of the components of 
the equations of motion, (12. 22 p . 

The map (ll.4j) cannot be transformed in a manner that puts into a canonical form without 
disrupting the deterministic normal form. In §3.3l and §3. 41 we showed that with certain reasonable 
assumptions, the shape of G has little effect on the resulting invariant density. 

Numerical simulations of (ll.4p revealed invariant densities that roughly conformed to one of 
four descriptions. If, in the absence of noise (e = 0) there exists an attracting period-n solution 
with no points on the switching manifold, then for sufficiently small e > the invariant density is 
approximately the sum of n Gaussians scaled by -, §4.11 The associated covariance matrices were 
computed in §3.21 In the same situation with larger values of e, iterates fall near the switching 
manifold with higher probability. Since the square-root singularity in (jl.4p induces a strongly 
nonlinear stretching of phase-space near the switching manifold, the Gaussian approximation 
worsens as e is increased, §4.21 With still larger e, the invariant density may be highly irregular 
but of a shape amenable to a simple explanation, §4.31 When attractors coexist in the absence 
of noise, for e > the invariant density is comprised of pieces about each deterministic attractor, 
§4.41 The rough magnitude of e for which the effect of the square-root singularity becomes 
significant, is influenced by the strength of attraction of the deterministic period- n solution and 
the proximity of the solution to the switching manifold. 

Given the above descriptions of the numerically computed invariant density at various fixed 
values of /i, we can now explain the dominate features of Figs. |2] and [3j We note that transi- 
tions in Figs. |2] and |3] are P -bifurcations [65j at which the invariant density changes character 
fundamentally. The sequences of transitions have not been described before as the underlying 
deterministic system is piecewise-smooth. 

In Fig. 121 for fi < —0.001, say, the invariant density is an approximate Gaussian centred at 
the fixed point, (II. 3p . Note that the width of the invariant density about the non-impacting 
solution is practically constant leading up to the grazing bifurcation, which is not the typical 
case for smooth maps [SJ EB]. For fi < 0, with \fi\ < O(e), the fixed point is close enough to the 
switching manifold that iterates repeatedly lie in the right half-plane. Due to the square-root 
singularity, the invariant density is non- Gaussian and 0(-\fe) in size. 

For small /i > 0, there is an irregular invariant density of size 0(^fe)\ roughly periodic 
behaviour at high periods is not seen. Over a wide interval of intermediate values of /j, the 
invariant density is well approximated by four Gaussians. As fi increases across this interval, 
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the approximation improves because the strength of attraction of the deterministic period-4 
solution increases with [i. The four-piece invariant density transforms to a widely spread density 
as the deterministic period-4 solution collides with the switching manifold in a border-collision 
bifurcation at /J, « 0.014. In the absence of noise, just beyond this bifurcation, there is a chaotic 
attractor. With noise the invariant density has a complicated wide-spread structure. For larger 
/j, the invariant density is well-approximated by three Gaussians. 

Fig. [3] displays much the same behaviour, but in addition iterates traverse neighbourhoods 
of coexisting deterministic attractors. Note that in some cases iterates only appear to surround 
a single attractor (e.g. for \i = 0.2). This occurs when the attractor is strongly attracting and 
its basin of attraction is large, relative to the noise amplitude, and the number of iterates of the 
numerical simulation is relatively small. 

It remains to investigate grazing bifurcations of vibro-impacting systems by studying a model 
with noise applied only to impact events. The formulation (11 .4p does not allow for an investigation 
into rapidly occurring, noise-induced multiple impacts. Many other bifurcation scenarios in 
piecewise-smooth systems await a detailed investigation in the presence of noise, such as sliding 
bifurcations. 



A The stochastic Nordmark map for the linear oscillator 

Here we derive the parameters, r, 5, \i an d @, of the stochastic Nordmark map, (II .4p . for the 
impact oscillator described in §2.31 We first let 
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(t — tgraz) mod 2"7T 



F — F, 



graz j 



t 



graz 



tan 



k 



(A.1) 



(A.2) 



is the time in (f , 7r) at which the impact occurs for the grazing periodic orbit, so that the system 
takes the general form (12.91) and satisfies the related assumptions (except phase space is 1R 2 x S 1 
instead of M 3 , but this is inconsequential). From the deterministic solution to (12. 23 p . we find 
that the global map, G, is given by (12. 6p with 
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in the discontinuity map, D, (12.51) . From the coordinate change (12. 7p . we obtain 



(A.3) 
(A.4) 
(A.5) 



r = 2e 2na cos(2vr/3) , 
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X = sgn(ai 2 c) 



(A.6) 
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where a ± i/3 are the eigenvalues of J, and specifically, 



P 



b 2 

osc 

4 



(A.7) 



It remains to use (12. 22ft to compute the covariance matrix, G, of the noise, £. Here we evaluate 
all terms of Q222), for QZ23) . In the form (EH]), we have 



/ (L) (u;^) 



fl(u; 7/) 



-k osc (u + 1) - 6 OSC V + fa + Fgraz) COs(u> + t graz ) 
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(A.8) 
(A.9) 



so by (E21 and ([Aj), ([0, 0, 0] T ; 0) = [0, - 7l , Cl] t = [0, -1, 1] T . Then (EHD yields 



e(t s) 
e 2 J 2 (t " s) 




(A.10) 



where e^- s ^ denotes the (i, j)-element of e J ^ S K Finally, the evaluation of (I2.16P at = 2ix, 
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using (IA.4P and (1A.10D . gives 
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and wi 3 = u 2 3 = uj 33 = 0. 
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